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1  Introduction 

In  this  lecture  we  shall  focus  on  the  dynamics  of  individual  phytoplankton  cells  and  on 
simple  models  to  describe  their  motion  in  the  fluid  flow.  Phytoplankton  cells  are  small, 
with  dimensions  up  to  a  few  hundred  microns.  Individual  phytoplankton  cells  immersed  in 
water  are  advected  by  the  surrounding  flow,  even  though  some  species  can  resist  passive 
advection  by  various  means,  such  as  varying  their  buoyancy  and/or  swimming.  In  addition, 
most  plankton  cells  are  slightly  heavier  than  water  and  are  thus  subject  to  gravitational 
sinking  (i.e.,  they  are  negatively  buoyant). 

If  we  want  to  describe  the  dynamics  of  a  phytoplankton  cell  advected  by  a  flow,  we  need 
to  write  the  equations  of  motion  for  the  cell.  However,  real  phytoplankton  cells  usually 
have  a  very  complex  shape  (which  can  change  with  time),  are  non-homogeneous  in  their 
interior,  can  vary  their  buoyancy  and  can  emit  substances  that  change  the  local  fluid  vis¬ 
cosity.  In  addition,  they  swim,  grow,  reproduce  (split),  and  eventually  die  or  are  eaten  by 
zooplankton  (small,  usually  multicellular  organisms  which  feed  on  algae).  For  fluffy  bodies 
with  spines,  tails  and  flagella  and  an  active  life  such  as  phytoplankters,  the  equations  of 
motion  are  rather  complicated,  and  one  should  resolve  the  Navier-Stokes  equations  with 
the  correct  boundary  conditions  on  the  mobile  surface  of  the  little  plankter. 

In  this  lecture  we  follow  a  traditional  path  for  the  GFD  school,  and  simplify  -  even 
grossly  oversimplify  -  the  problem,  by  considering  the  dynamics  of  infinitesimally  small, 
spherical  heavy  particles,  with  fixed  shape  and  density,  which  are  passively  advected  by  a 
surrounding  fluid  flow,  and  we  call  them  “impurities.”  Of  course,  our  impurities  are  but  a 
far  image  of  phytoplankton  cells  in  the  ocean. 

2  Plankton  and  turbulence:  a  matter  of  size 

Most  fluid  flows  in  the  ocean  and  in  lakes  are  turbulent,  especially  in  the  upper  mixed  layer 
where  phytoplankton  thrive.  Thus,  on  large  scales,  the  world  in  which  plankton  cells  live  is 
a  turbulent  environment.  However,  individual  cells  are  small,  usually  smaller  than  the  Kol¬ 
mogorov  scale  -  the  length  scale  which  marks  the  transition  from  laminar  flow  at  smaller 
scales  to  turbulent  dynamics  at  larger  scales.  For  this  reason,  individual  phytoplankters 
experience  a  viscous,  laminar  and  non-turbulent  flow  in  their  immediate  surroundings. 
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The  motion  of  sea  or  lake  water  under  “normal”  conditions  (ignoring  compressibility, 
thermodynamics,  phase  transitions  and  other  complications)  can  be  described  by  the  in¬ 
compressible  Navier-Stokes  equations, 

Du  „  _ tj  „ 

p— = -Vp  + pgz  +  F  (1) 

where  u(x,  t)  is  the  fluid  velocity,  x  is  the  spatial  coordinate,  t  is  time,  p  is  the  fluid  density, 
p  is  pressure,  p,  is  the  dynamic  viscosity,  g  is  the  gravitational  acceleration,  z  is  the  unit 
vector  in  the  z  direction,  F  is  a  generic  body  force  (if  present)  and  the  substantial  (or 
Lagrangian)  derivative  can  be  expressed  as  Du/ Dt  =  du/dt  +  u  •  Vu.  We  have  fixed  the  z 
coordinate  to  point  vertically  upwards.  If  p  is  assumed  to  be  constant,  the  gravitation  term 
can  be  absorbed  in  the  definition  of  pressure. 

Equation  (1)  is  complemented  by  the  incompressible  version  of  the  continuity  equation, 

V  •  u  =  0  ,  (2) 

and  in  this  way  we  get  a  complete  description  of  the  fluid  flow,  four  equations  in  four  vari¬ 
ables.  Of  course,  to  (try  to)  solve  the  equations  we  also  need  boundary  and  initial  conditions. 

The  Reynolds  number  can  be  defined  as  the  ratio  between  the  strength  of  the  nonlinear 
term,  p\u-  Vu|,  and  that  of  the  viscous  term,  |^V2u|.  If  we  call  U  and  L  the  typical  scale  of 
velocity  and  length,  and  scale  time  by  L/U,  we  get  that  the  Reynolds  number  is  R e=UL/u 
where  v  =  p,/ p  is  the  kinematic  viscosity.  When  the  Reynolds  number  is  smaller  than  one, 
i.e. ,  the  velocity  and/or  the  length  scale  are  sufficiently  small,  the  flow  is  laminar  and  viscos¬ 
ity  dominates.  For  Re-^  0,  we  can  discard  the  inertial  term  in  the  Navier-Stokes  equations 
and  find  viscosity-dominated  special  solutions.  When  ReS>  1,  the  flow  becomes  turbulent 
and  the  fluid  velocity  varies  intermittently  and  erratically  in  time  and  space,  spanning  a 
wide  range  of  values. 

The  simplified  form  of  the  Navier  Stokes  equations  given  by  eq.(l,2)  with  constant 
density  describes  the  cleanest  case  of  turbulence,  i.e.,  homogeneous  and  isotropic  three- 
dimensional  turbulence.  Physically,  this  type  of  flow  is  a  good  model  for  turbulence  for 
scales  up  to  a  few  meters,  a  reasonably-sized  portion  of  the  world  for  a  phytoplankton  cell. 

To  quantitatively  describe  this  type  of  turbulence,  we  assume  that  the  energy  input  in 
the  flow  takes  place  at  a  certain  macro  length  scale,  Iq.  and  dissipation  takes  place  mainly 
at  a  very  small  scale,  say  Id-  In  the  scale  range  between  Iq  and  Id ,  energy  gets  transferred 
to  smaller  and  smaller  length  scales  by  the  breakdown  of  the  larger  eddies  into  smaller  ones. 
Here,  an  “eddy”  is  a  generic,  usually  short-lived,  vortical  structure  in  the  three-dimensional 
turbulent  flow,  and  it  is  characterized  by  a  given  (approximate)  length  scale.  If  Iq  3>  Id 
and  no  forcing  is  present  between  these  two  scales,  then  we  see  the  emergence  of  an  inertial 
range  where  the  dynamics  is  independent  of  dissipation  and  forcing.  Assuming  that  the 
spectral  energy  transfer  rate  is  constant  in  the  inertial  range,  so  to  obtain  a  statistically 
stationary  energy  spectrum,  and  using  simple  scaling  arguments,  one  can  derive  the  form  of 
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the  energy  spectrum  in  the  inertial  range,  E(k),  which  obeys  the  celebrated  Kolmogorov’s 
-5/3  law, 

E(k)  ~  e2/3A;-5/3  (3) 

where  k  is  the  wavenumber  (i.e.,  2tt/\  where  A  is  the  wavelength)  and  e  is  the  energy  dis¬ 
sipation  rate  per  unit  volume. 


Figure  1:  Energy  spectrum  of  a  three-dimensional,  homogeneous  and  isotropic  turbulent 
flow. 

The  length  scale  at  which  diffusion  becomes  important  is  called  the  Kolmogorov  scale, 
lv ,  defined  as 

lu  ~  k-1  ~  (vZ/e)1'4  .  (4) 

At  scales  smaller  than  lu,  the  flow  is  laminar  (but  possibly  unsteady).  In  the  ocean,  lv 
ranges  between  6  x  10-5  m  and  0.01  m.  Since  most  phytoplankton  cells  are  smaller  than  a 
few  tens  of  microns  in  size,  they  do  not  experience  a  turbulent  flow  but  are  rather  exposed 
to  a  viscous  laminar  flow. 


The  dominance  of  viscosity  makes  nutrient  uptake  by  phytoplankton  to  be  severely 
limited  by  molecular  diffusion.  For  a  plankton  cell,  nutrient  delivery  happens  in  two  stages 
-  nutrient  delivery  to  the  cell  membrane  and  transport  across  the  membrane.  When  the 
latter  process  is  faster  than  the  first  one  (as  often  happens),  the  nutrient  concentration  at 
the  cell  surface  drops  significantly  below  its  far-held  value.  For  a  molecular  diffusivity  D 
and  a  plankton  cell  of  radius  a,  the  concentration  of  nutrient  in  the  proximity  of  the  cell 
wall  is  governed  by  the  equation 


=  /;v2.V 

dt 


(5) 


whose  spherically  symmetric,  steady  state  solution  is 

N(r)  =  -(N0-Noo)  +  Noo  (6) 

r 

where  Nq  is  the  concentration  at  the  cell  surface  and  is  the  far-held  concentration.  The 
diffusional  hux  to  a  spherical  cell  is  then  computed  as 


Qd  =  47rDa(lV00  -  Nq)  . 


(7) 
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When  the  uptake  capacity  of  the  cell  exceeds  the  rate  of  nutrient  supply  by  diffusion,  a 
nutrient-depleted  layer  is  formed  near  the  cell.  The  time  taken  by  diffusion  for  transporting 
the  nutrient  across  this  layer  can  be  estimated  as 


Td  = 


D 


(8) 


In  water  at  standard  upper  ocean  temperatures,  for  innnotile  cells  with  a  <  lO^um  the  nu¬ 
trient  limitation  can  be  severe. 


Given  that  diffusion  is  not  sufficient  to  bring  enough  nutrients  to  the  cell  surface,  es¬ 
pecially  for  smaller  organisms,  phytoplankters  have  found  ways  to  overcome  diffusion  lim¬ 
itation.  Most  of  these  strategies  are  based  on  the  ability  of  phytoplankters  to  move  along 
trajectories  which  are  not  fluid  particle  paths,  thus  allowing  them  to  explore  new  regions  of 
the  fluid  which  are  not  nutrient  depleted  yet.  The  differences  between  cell  trajectories  and 
fluid  parcel  trajectories  can  be  generated,  for  example,  by 

1.  Swimming  (as  many  flagellates  do); 

2.  Sinking  (as  diatoms  do); 

3.  Exploiting  the  fact  that  plankton  cells  have  different  density  from  the  fluid  and  have 
finite  size. 

In  the  rest  of  this  Lecture,  we  shall  focus  on  the  last  mechanism  and  study  how  heavy 
impurities  with  finite  size  (our  gross  model  for  plankters)  are  advected  by  a  fluid  flow. 


3  Lagrangian  advection 

The  motion  of  an  advected  particle  with  the  same  density  of  the  fluid  and  vanishingly  small 
dimensions  (i.e.,  a  fluid  parcel)  is  described  by  the  standard  Lagrangian  equations 

dX 

~dt  ~  U(X’  *)  <9> 

where  X  is  the  particle  position  at  time  t  and  u(X,  t )  is  the  fluid  velocity  at  the  particle 
position.  Notice  that  equation  (9)  is  non-Newtonian,  in  the  sense  that  it  does  not  obey 
the  standard  Newton’s  second  law  of  motion.  The  reason  is  that  the  particle  has  no  inertia 
with  respect  to  the  fluid,  and  thus  we  equate  the  Lagrangian  particle  velocity  dX/dt  to  the 
Eulerian  velocity  u(X,f). 

Impurities  with  finite  size  and  finite  inertia  do  not  behave  like  fluid  particles.  An  impu¬ 
rity  with  density  pp  moving  with  velocity  V  =  dX /dt  in  a  flow  field  u  of  a  fluid  with  density 
Pf  experiences  forces  due  to  gravity,  F9,  Stokes  drag,  Fs,  and  the  effect  of  the  surrounding 
fluid  (mainly  pressure),  F f,  expressed  as 

1  9  upf  Du 

Fg  =  -g{pP~  pf)z,  Fs  =  -j— -^pfua(\-u)  =  — 2^-(V-u),  Ff  =  pf— 

(10) 
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where,  as  before,  a  is  the  radius  of  the  impurity. 


Thus,  using  Newton’s  second  law  equating  force  to  mass  times  acceleration,  we  can  write 
for  the  motion  of  the  impurity 

dV 

=  Fj  +  Fs  +  Fg  .  (11) 

Rewriting  equation  (11)  we  have 

W  =  »)-('-%*  (12) 

where  5  =  Pf/Pp  and  ra  =  2a2/  (95V).  The  time  scale  ra  measures  the  relaxation  time  of 
the  particle  velocity  to  the  local  fluid  velocity  and  6  determines  the  nature  of  the  buoyant 
particle  (6  <  1  and  5  >  1  indicate  particles  which  are  heavier  or  lighter  than  the  fluid 
respectively). 


In  still  fluid,  where  u 
impurity  as 


0,  one  can  define  the  (downward)  terminal  velocity  of  the  heavy 


2  Pp _ Pf_ 

PfV 


(13) 


For  example,  an  impurity  with  5  =  0.97  and  a  =  30/rm,  as  appropriate  for  a  plankton  cell, 
the  terminal  velocity  W  is  a  few  meters  per  day. 


It  is  important  to  stress  the  difference  between  the  two  time  derivatives  present  in  equa¬ 
tion  (14).  While  D/Dt  =  d/dt  +  u  •  V  measures  the  rate  of  change  along  the  path  that  a 
fluid  particle  would  follow  locally,  d/dt  =  d/dt  +  V  •  V  denotes  the  time  derivative  while 
tagging  along  the  heavy  impurity. 


/ 

/ 

I 


Streamline 

I 

I 
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Figure  2:  A  graphical  representation  of  the  meanings  of  d/dt  and  D/Dt 


Before  proceeding  further,  we  make  eq.(12)  non  dimensional  and  write 

d' V  Du  1  ~  x\  1  ~ 

dt  Dt  SV  J  v  '  Fr2 


(14) 
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where  the  tilde  is  used  to  denote  dimensionless  variables.  In  the  non-dimensional  equation, 
St  =  TaU/L  and  Fr  =  U/y/gL  are  the  Stokes  and  Froude  numbers  respectively  and  L 
and  U  are  typical  length  and  velocity  scales.  The  Stokes  number  St  represents  the  ratio 
between  the  particle  relaxation  time  scale,  r0,  and  the  dynamical  time  scale  L/U,  whereas 
Fr  quantifies  the  ratio  of  the  buoyancy  time  scale,  l/^/g/L,  to  the  dynamical  time  scale.  In 
the  following  we  will  drop  the  tilde  and  denote  the  non-dimensional  variables  by  the  same 
notation  as  the  dimensional  ones,  unless  a  distinction  needs  to  be  made. 


4  Impurity  advection  in  2D  stationary  flows 

In  this  section  we  study  the  dynamics  of  an  ensemble  of  passive  impurities  advected  by  a 
stationary,  divergence-free,  two-dimensional  (2D)  cellular  velocity  field  without  the  effects 
of  gravity  (we  can  think  of  a  flow  field  on  a  horizontal  plane).  This  example  allows  for 
stressing  the  differences  between  the  dynamics  of  fluid  parcels  and  that  of  heavy  and  light 
impurities;  note  that  in  this  simple  case  one  can  rescale  the  Stokes  number  by  rescaling  time. 

We  first  consider  the  case  of  an  isotropic  streamfunction  =  2(cosx  +  cos  y).  The 
fluid  velocity  is  u  =  (u,v)  and  the  velocity  components  are  given  by  u  =  —dip/dy  and 
v  =  difr/dx.  In  a  2D  flow,  fluid  parcels  obey  the  equations 

dX  _  di[j  dY  _  ch/>  ,  . 

~dt  =  ~dY  ’  ~dt  =  OX  (  ’ 

which  is  formally  an  Hamiltonian  system  with  phase  space  (x,  y )  where  x  and  y  are  canoni¬ 
cally  conjugate  variables  and  the  stream  function  plays  the  role  of  the  Hamiltonian.  In  this 
case,  one  has  conservation  of  phase-space  area  and  no  preferential  fluid  particle  concentra¬ 
tion  is  possible.  For  a  stationary  flow,  fluid  parcels  follow  the  streamlines;  since  the  system 
is  two-dimensional,  there  is  no  chaotic  dynamics  and  the  fluid  parcels  undergo  only  regular 
trajectories. 

The  situation  is  different  for  passive  impurities,  for  which  the  phase  space  (x,y,Vx,Vy) 
is  four-dimensional  and  the  dynamics  is  dissipative.  In  the  case  of  particles  lighter  than  the 
fluid,  6  >  1,  linear  stability  calculations  show  that  the  fixed  points  in  (xn,  ym)  =  (2n7r,  2rmr) 
and  {x'n.  y'm)  =  ((2ra  +  l)7r,  (2ra  +  l)7r)  ( n,m  being  integers)  are  stable  and  attracting,  while 
the  fixed  points  (. xn,ym )  =  (2n7r,  (2m  +  l)7r)  and  (. x'n,y'm )  =  ((2 n  +  l)7r,2m7r)  are  sad¬ 
dles.  For  heavy  particles  (6  <  1),  all  these  fixed  points  are  unstable.  In  this  case,  the 
impurities  undergo  chaotic  trajectories,  hopping  from  one  flow  cell  to  another,  and  display 
non-turbulent  diffusion,  even  in  the  case  of  a  stationary  velocity  field  [3]. 

Figure  4  shows  how  an  initially  homogeneous  distribution  of  heavy  and  light  impurities 
evolves  in  the  prescribed  flow.  The  short-time  dynamics  is  very  similar  for  both  light  and 
heavy  particles.  An  interesting  feature  is  the  transformation  of  an  initially  uniform  particle 
distribution  into  a  distribution  characterized  by  sharp-edged  caustics;  note  that  caustics 
cannot  form  for  fluid  particles,  owing  to  the  conservation  of  phase-plane  area  associated 
with  the  Hamiltonian  nature  of  the  system.  This  constraint  does  not  act  for  impurites, 
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Figure  3:  The  stationary  two-dimensional  streamfunction  ip  =  2(cosx  +  cos  y)  and  the 
associated  velocity  field. 

owing  to  the  dissipative  nature  of  the  equations  describing  impurity  motion.  For  light 
particles,  the  caustics  spiral  around  the  stable  fixed  points.  Hence,  for  a  stationary  two 
dimensional  flows,  light  particles  tend  to  concentrate  at  the  stationary  points  of  the  flow; 
for  large  enough  times  the  particles  are  practically  absent  in  the  remaining  regions  of  the 
flow.  On  the  contrary,  heavy  particles  exhibit  a  very  different  behavior:  After  the  initial 
formation  of  caustics,  the  heavy  particles  undergo  diffusion  along  the  separatrices  of  the  ad¬ 
verting  flow,  owing  to  the  chaotic  nature  of  the  heavy  particle  dynamics  for  the  parameter 
values  adopted  here. 

Figure  5(a)  shows  the  trajectory  of  an  individual  heavy  impurity  in  the  cellular  flow 
field.  The  particle  alternates  between  spiralling  inward  around  one  of  the  (now  unstable) 
fixed  points  and  escaping  along  one  of  the  unstable  directions,  with  a  quasiballistic  mo¬ 
tion.  Although  this  alternance  of  trapping  and  quasiballistic  motion  could  be  reminiscent 
of  a  Levy  flight,  for  the  system  under  study  the  long-time  probability  distribution  of  the 
particle  velocity  is  Gaussian,  and  it  is  consistent  with  a  standard  Brownian  diffusion  process. 

We  can  compute  the  diffusion  coefficient  for  an  ensemble  of  N  advected  impurities  as  a 
function  of  e  =  1  —  5, 


Dx  =  9^V  X  t Xi ®  ~  Xi(°)]2  (16) 

2—1 

1  N 

Dy  =  &  9 Jn  X  °)]2  •  (17) 

2—1 

Figure  5(b)  shows  the  dependence  of  Dx,Dy  on  e.  Since  the  flow  is  isotropic,  the  two  dif¬ 
fusion  coefficients  are  equal.  In  general,  no  anomalous  (non-Brownian)  diffusion  has  been 
observed  for  this  system. 


Ill 


In  the  case  of  a  weakly  anisotropic  velocity  field,  one  observes  some  new  effects.  In 
the  following,  we  consider  the  flow  described  by  the  streamfunction  ip  =  A  cos  x  +  B  cos  y, 
where  A  ^  B.  If  one  looks  at  the  streamfunction  (figure  6),  it  is  evident  that  there  are 
open  “channels”  (with  our  choice  of  parameters,  in  the  horizontal  direction)  where  the  fluid 
particles  undergo  ballistic  motion.  These  are  analogous  to  the  libration  orbits  in  the  phase 
space  of  a  nonlinear  pendulum  (the  closed  orbits  would  correspond  to  periodic  oscillations). 
Thus,  fluid  particles  are  either  trapped  inside  the  cellular  regions  or  undergo  ondulating 
ballistic  motion  in  the  channels. 

A  different  story  is  observed  for  advected  impurities.  As  for  the  isotropic  flow,  light 
particles  sooner  or  later  get  trapped  at  the  stagnation  points.  On  the  other  hand,  heavy 
impurities  spend  a  long  time  in  quasiballistic,  motion  in  the  channels,  and  temporary  trap¬ 
ping  in  the  cellular  regions  is  followed  by  jumping  off  to  another  channel.  Figure  7  shows 
that  an  intially  uniform  particle  distribution  tends  to  become  spread  along  the  channels, 
with  strong  anisotropy. 

In  this  flow,  the  diffusion  coefficients  along  the  x  and  y  directions  are  very  different.  Fig¬ 
ure  8(b)  reveals  that  Dx ,  the  diffusion  coefficient  along  the  channels,  diverges  with  e  — >  0, 
Dx  ~  e-°-6  while  Dy  ~  e1-5.  Thus  a  small  anisotropy  in  the  Eulerian  flow  {A  =  2,  B  =  2.2) 
leads  to  diffusion  coefficients  that  differ  by  several  orders  in  magnitude  in  the  case  of  par¬ 
ticles  which  are  only  slightly  heavier  than  the  fluid  {5  =  0.9).  This  could  have  interesting 
consequences  for  the  motion  of  plankton  cells  in  weakly  anisotropic  flows,  as  the  small  den¬ 
sity  difference  between  the  plankton  cells  and  the  water  could  result  in  a  strongly  anisotropic 
diffusion  of  the  particles. 


5  Neutrally  buoyant  impurities  with  finite  size 

Often,  advected  impurities  have  a  density  that  is  rather  close  to  that  of  the  fluid  (as  hap¬ 
pens,  for  example,  with  most  plankton  cells),  but  have  a  finite  size,  i.e.,  they  are  not  fluid 
elements.  To  explore  the  dynamics  of  these  particles,  we  consider  the  case  of  neutrally 
buoyant  particles  with  5  =  pf/pp  =  1  and  finite,  albeit  small,  radius  a.  The  equation  for 
particle  motion  in  this  case  becomes: 


where  now 


(TV 

dt 


Du 

~Dt 


St 


2 

9 


(18) 


If  we  make  the  approximation  that  the  time  derivative  along  the  trajectory  of  an  impurity 
and  that  along  the  trajectory  of  a  fluid  parcel  are  the  same,  i.e.  d/dt  =  D /Dt,  the  problem 
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becomes  very  simple.  In  this  case,  equation  (18)  can  be  written  as 

|,V-u)  =  -i(V-u,. 

Hence 

V  =  u  +  (V  —  u)0  exp^—  ^tj  , 

where  (V  —  u)0  denotes  the  initial  value  of  (V  —  u).  From  this  we  could  infer  that  even 
if  we  release  a  particle  with  an  initial  velocity  which  differs  from  that  of  the  fluid,  after 
a  transient  the  particle  velocity  will  tend  to  the  fluid  velocity,  implying  that  a  neutrally 
buoyant  particle  will  behave,  after  an  initial  transient,  as  an  ideal  tracer. 


In  reality,  the  above  conclusion  is  not  correct,  as  the  problem  is  complicated  by  the  fact 
that  d/dt  and  D/Dt  are  not  the  same.  Equation  (18)  should  then  be  written  as 

^  (V  -  u)  =  -  (V  -  u)  ■  Vu  -  -i  (V  -  u) .  (19) 

Defining  M  =  V  —  u,  equation  (19)  can  be  rewritten  as 


dM 

dt 


-M-  Vu-  —  M 

St 


-|J  +  ^ 


•M, 


(20) 


where  I  is  the  identity  matrix  and  J  is  the  velocity  Jacobian, 


f dxu  dyu\ 
\dxv  dyv]  ’ 


u(x) 


If  we  diagonalise  the  matrix  we  obtain 


dMd 

dt 


■Md. 


From  this  expression,  one  can  obtain  relevant  information  on  the  behavior  of  the  advected 
impurities.  For  a  two-dimensional  stationary  flow,  the  difference  between  the  fluid  velocity 
and  the  particle  velocity  evolves  in  time  depending  on  the  value  of  ±A  —  1/St.  In  particular, 
the  particle  velocity  can  diverge  from  the  fluid  velocity  when 


from  which  one  sees  that  the  stability  of  the  particle  trajectories  (with  respect  to  the  fluid 
parcel  trajectories)  is  determined  by  the  balance  between  the  eigenvalues  of  the  fluid  ve¬ 
locity  field  and  the  Stokes  number.  When  the  local  value  of  A2  (a  property  of  the  Eulerian 
advecting  flow)  is  sufficiently  large,  the  neutrally  buoyant  impurity  can  follow  a  trajectory 
that  is  different  from  that  of  a  fluid  parcel  placed  at  that  point. 
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We  can  get  a  further  advancement  by  resorting  to  the  Okubo- Weiss  parameter  [7,  11] 

Q  =  -det  J  =  j  (S2  -  uj2)  ,  (21) 

where  S 2  =  ( dxu  —  dyv )2  +  (dyu  +  dxv )2  is  the  total  strain  (deformation)  and  u>  =  dv/dx  — 
du/dy  is  the  vorticity.  This  parameter  characterises  the  local  shear/rotational  properties 
of  the  flow  and  separates  the  fluid  domain  into  elliptic  and  hyperbolic  regions  [7,  11,  4]: 


Q 


<  0  rotation  —  dominated 
>  0  shear  —  dominated. 


Since  A  satisfies  det(J  —  AI)  =  0,  one  has  that 

A2  —  7YJ  +  det  J  =  0  .  (22) 

In  addition,  as  we  are  now  considering  an  incompressible  flow, 

V  •  u  =  dxu  +  dyv  =  Tr J  =  0.  (23) 

From  (22)  and  (23),  the  Okubo- Weiss  parameter  is  then  Q  =  A2.  If  Q  >  0,  A  is  real  and 
deformation  dominates,  as  around  hyperbolic  points,  whereas  if  Q  <  0,  A  is  imaginary  and 
rotation  dominates,  as  in  the  proximity  of  elliptic  points.  Thus,  we  see  that  hyperbolic 
regions  can  lead  to  a  divergence  between  the  fluid  particle  paths  and  the  trajectories  of  the 
impurities,  if  they  are  “hyperbolic  enough”  (i.e.,  the  value  of  A2  is  large  enough).  As  a 
result,  neutrally  buoyant  impurities  will  be  ejected  from  hyperbolic  regions  and  will  wander 
around  till  they  end  up  in  elliptic  regions,  were  they  will  follow  the  same  trajectory  of 
a  fluid  parcel.  This  potentially  leads  to  a  preferential  concentration  of  neutrally  buoyant 
impurities  in  elliptic  regions,  with  a  biassed  sampling  of  the  advecting  flow. 

5.1  Example  1:  Time  independent  flow 

To  illustrate  the  effects  of  St  and  Q  on  the  dynamics  of  a  neutrally  buoyant  particle,  we 
consider  the  simple  incompressible  two-dimensional  model  flow  defined  by  the  stream  func¬ 
tion  i/j(x,y)  =  A  cos  x  cos  y  [1]. 

Figure  9(a)  shows  the  contour  lines  of  the  stream  function  and  of  the  Okubo- Weiss 
parameter  Q.  Figure  9(b)  shows  the  trajectory  of  a  neutrally  buoyant  particle  whose  initial 
velocity  is  different  from  that  of  the  fluid  at  the  impurity  position.  The  difference  between 
the  velocity  of  the  impurity  and  the  fluid  velocity  is  amplified  in  highly  hyperbolic  regions, 
where  impurity  trajectories  and  fluid  parcel  trajectories  separate  exponentially.  Fig.9(c) 
shows  that  the  neutrally  buoyant  impurity  has  a  long  and  complicated  trajectory,  wandering 
between  different  cells  before  it  comes  across  the  vicinity  of  an  elliptic  point  and  settles  in 
one  of  the  low-Q  regions.  The  time  evolution  of  the  difference  between  the  impurity  velocity 
and  the  fluid  velocity  at  the  impurity  position  is  shown  in  Fig. 9(d).  Once  the  impurity  has 
settled  in  an  ellpitic  region,  the  velocity  difference  between  the  fluid  and  the  impurity  tends 
to  zero  an  no  further  instability  is  observed.  However,  in  their  wanderings  impurities  leave 
hyperbolic  regions,  leading  to  a  preferential  concentration  of  neutrally  buoyant  impurities 
in  the  elliptic  regions  of  the  flow. 
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5.2  Example  2:  Time  periodic  flow 

The  time-dependent  flow  defined  by  the  stream  function  tp(x,  y)  =  A  cos(.x  +  B  sin  ujt)  cos  y 
shows  interesting  features  [1] .  As  in  a  typical  non-stationary,  two-dimensional  non-integrable 
Hamiltonian  system,  there  are  regions  of  phase  space  where  fluid  parcels  undergo  chaotic 
motion.  Figure  10(a)  shows  the  Poincare  section  of  the  chaotic  trajectory  of  a  fluid  particle, 
here  obtained  by  plotting  the  parcel  position  at  times  t  =  0,  T,  2 T, ...  where  T  =  2i t/lo.  The 
trajectory  visits  a  large  part  of  the  domain  but  it  is  bounded  by  two  separatrices  parallel  to 
the  x  direction,  which  forbid  fluid  particle  motion  across  them  (remember  that  the  temporal 
variability  is  only  in  the  velocity  component  along  x). 

The  trajectory  of  a  neutrally  buoyant  impurity,  shown  in  Fig.  10(b),  indicates  that  the 
impurity  can  wander  on  a  larger  portion  of  the  domain,  but  in  the  end  it  converges  to  a 
non-chaotic  trajectory  around  an  elliptic  point.  From  that  moment,  the  impurity  behaves  as 
a  fluid  particle.  However,  the  impurity  has  selected  a  specific  sub-domain  in  space  without 
ever  coming  back  close  to  its  initial  conditions,  at  variance  with  what  would  have  happened 
for  a  Hamiltonian  system.  The  equations  of  motion  for  a  neutrally  buoyant  impurity  are 
dissipative,  and  there  is  no  conservation  of  phase-space  volume.  For  completeness,  Fig.  10(c) 
shows  the  difference  between  the  velocity  of  the  impurity  and  that  of  the  fluid,  as  a  function 
of  time. 

5.3  An  application  to  two-dimensional  turbulence 

The  results  discussed  above  have  interesting  implications  for  two-dimensional  turbulent 
flows  [1,  9].  We  integrate  the  two-dimensional  vorticity  equation  for  an  incompressible  fluid 
with  bi-harmonic  (hyper)viscosity, 

did 

—  +  uj)  =  -vhVAui,  (24) 

where  i/j  is  the  stream  function,  oj  =  V2?/’  is  the  vorticity,  u  =  —dyip  and  v  =  dxip  are  the 
velocity  components  and  the  hyperviscosity  parameter  is  t'/,  =  2.5  x  10~7.  To  integrate  the 
equation,  we  use  a  pseudo-spectral  method  with  5122  collocation  points  in  a  doubly  periodic 
domain;  here  we  consider  a  freely-decaying  case  with  no  external  forcing.  In  this  flow,  an 
initially  homogeneous  distribution  of  neutrally  buoyant  impurities  evolves  in  time  towards 
a  distribution  where  the  impurities  are  concentrated  in  the  inner  part  of  vortices,  where 
Q  <  0,  as  shown  in  Fig.  11.  This  may  have  interesting  consequences  when  we  design  fluid 
experiments  using  neutrally  buoyant  tracer  particles,  as  one  should  verify  that  the  tracers 
really  provide  an  unbiassed  sampling  of  the  velocity  field. 

6  The  role  of  rotation 

The  presence  of  a  rotation  introduces  new  effects  in  the  dynamics  of  the  heavy  impurities 
[9].  When  the  (apparent)  Coriolis  and  centrifugal  forces  are  acting  on  an  impurity,  the 
dimensional  equation  for  the  impurity  motion  becomes 

— —  =  5— - (V  —  u)  —  (1  —  6)  gz  —  2f2  x  (V  —  <5u)  +  |St|2  R  (1  —  S) ,  (25) 

at  Dt  Ta 
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where  6  and  ra  are  defined  as  before  and  R  is  the  distance  of  the  particle  from  the  axis 
of  rotation.  Notice  that  the  term  —2 SCI  x  u  had  to  be  subtracted  from  the  Coriolis  term 
acting  on  the  impurity  because  it  was  already  taken  into  account  in  the  expression  dDu/Dt. 
Thus,  the  Coriolis  term  acts  on  the  density-weighted  difference  between  the  velocity  of  the 
impurity  and  that  of  the  fluid.  Non-dimensionalising  the  above  equation,  we  obtain 


dV 

dt 


(26) 


where  St  is  defined  as  before  and 


Ro  = 


U 

2RL 


is  the  Rossby  number.  Fig.  12  shows  a  distribution  of  8000  heavy  impurities  which  were  ini¬ 
tially  homogeneously  distributed  in  a  statistically  stationary  numerical  simulation  of  forced 
and  dissipated  2D  turbulence  [9].  The  numerical  integration  has  been  performed  in  the 
doubly  periodic  domain  [0,  27r]  x  [0,  27t]  with  128  x  128  collocation  points.  The  particles 
have  5  =  0.8,  St  =  0.01  and  Ro  =  1/320.  In  the  integration,  we  have  discarded  the  effect  of 
the  centrifugal  force.  The  heavy  impurities  concentrate  inside  the  centers  of  the  anticyclonic 
vortices,  leaving  both  the  turbulent  background  and  the  cores  of  cyclonic  vortices.  This  is 
due  to  the  fact  that  the  Coliolis  force  adds  up  to  the  inertial  force  to  eject  heavy  impurities 
from  cyclones,  whereas  it  works  against  the  inertia  to  push  impurities  towards  the  center  of 
anticyclones.  When  Ro  is  small  enough,  the  Coriolis  force  overcomes  the  centrifugal  effect 
of  inertia  and  the  impurities  concentrate  in  the  cores  of  anticyclonic  vortices.  When  the 
vortices  are  especially  intense,  they  generate  a  strong  centrifugal  force  due  to  the  particle  in¬ 
ertia.  The  strength  of  the  centrifugal  push  by  the  vortex  becomes  weaker  when  the  distance 
from  the  vortex  center  increases.  In  this  case,  the  Coriolis  force  balances  inertia  at  some 
distance  from  the  vortex  center  and  the  impurities  concentrate  in  an  annular  distribution 
inside  the  vortex. 


As  a  final  example,  we  momentarily  leave  the  world  of  plankton  and  discuss  the  (again 
idealized!)  motion  of  dust  grains  in  the  solar  nebula  [2].  In  this  case,  <5  ~  10-8  so  that  we 
may  may  discard  the  term  5  D\ /Dt.  The  impurity  equation  is  then  reduced  to 

^  (V  -u)  -  V<F  -  2fi  x  V  +  |S1|2R,  (27) 

dt  Ta 

where  $  is  a  gravitational  potential.  In  the  specific  case  of  particles  spinning  around  the 
protosun  in  the  solar  nebula  (a  rarefied  fluid  composed  mainly  by  hydrogen)  the  equation 
of  motion  of  a  dust  particle  becomes 

—  (V-u)  -2U  x  V-  ^R+  |II|2R.  (28) 

dt  te  R 


Here,  G  and  M  are  the  gravitational  constant  and  the  mass  of  the  protosun,  respectively. 
For  the  rarefied  conditions  of  the  protoplanetary  nebula,  the  friction  time  scale  te  is  set  to 
have  the  form  te  oc  a/pf,  where  a  is  the  radius  of  the  particle.  This  expression  of  the  drag 
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coefficient  corresponds  to  what  is  called  the  Epstein  regime  for  a  rarefied  gas. 

To  proceed  with  the  study,  we  consider  a  simplified  representation  of  the  protosolar  neb¬ 
ula  by  using  a  two-dimensional  turbulent  flow  subject  to  a  circular ly-symmetric  Keplerian 
shear,  physically  generated  by  the  presence  of  a  central  massive  object  (the  protosun)  [2]. 
Figure  13(a)  shows  the  vorticity  field  at  sufficiently  large  time,  obtained  by  numerical  sim¬ 
ulation  of  this  system.  By  perturbing  the  initial  Keplerian  shear  (which  is  linearly  stable  in 
the  absence  of  magnetic  effects),  both  cyclonic  and  anticyclonic  vortical  perturbations  tend 
to  form.  However,  cyclonic  vortices  are  rapidly  destroyed  by  the  adverse  shear  -  or  do  not 
form  at  all  -  and  only  anticyclonic  vortices  survive  in  the  Keplerian  shear,  as  illustrated  in 
Figurel3(a).  As  already  seen  in  the  case  of  standard  two-dimensional  turbulence  discussed 
above,  the  heavy  impurities  tend  to  concentrate  inside  the  anticyclonic  vortices  as  time 
progresses.  Fig.  13(b)  shows  the  distribution  of  an  ensemble  of  heavy  impurities,  initially 
released  with  a  uniform  distribution,  in  the  flow  field  of  Fig.  13(a).  From  these  results,  one 
could  speculate  that  dust  grains  in  the  solar  nebula  could  be  concentrated  by  the  dynamics 
in  the  cores  of  anticyclonic  vortices,  where  the  subsequent  evolution  (collisions,  gravitational 
attraction)  could  more  easily  ensue. 


7  Conclusion 

The  dynamics  of  passively  advected  impurities  is  quite  different  from  that  of  fluid  parcels. 
From  a  mathematical  point  of  view,  the  phase  space  of  the  impurity  motion  has  twice  the 
number  of  dimensions  of  that  of  the  fluid  parcels,  since  it  is  described  by  particle  posi¬ 
tions  and  velocities  instead  of  just  particle  positions,  as  in  the  case  of  fluid  elements.  In 
addition,  the  presence  of  Stokes  drag  makes  the  system  dissipative,  with  the  possibility  of 
non-homogeneous  particle  distributions  in  space  and  preferential  particle  concentrations  in 
specific  regions  of  the  flow.  Physically,  this  means  that  advected  impurities  may  provide  a 
non-homogeneous  sampling  of  the  fluid  domain.  Even  neutrally  buoyant  impurities,  when 
their  size  is  large  enough,  can  behave  very  differently  from  fluid  parcels.  All  these  results 
indicate  that  small  objects  immersed  in  a  fluid  -  such  as  plankton  cells  in  water  or  dust 
grains  in  the  solar  nebula  -  can  undergo  complicated  dynamics  which  cannot  be  understood 
only  from  considering  the  motion  of  fluid  particles. 

Before  closing  this  lecture,  it  is  important  to  recall  that  so  far  we  have  simplified  the 
description,  neglecting  for  example  the  effect  of  unsteady  forces  associated  with  the  acceler¬ 
ation  of  the  impurities,  such  as  the  added  mass  term  [6]  and  the  Basset  force  [10].  These  two 
terms  become  important  especially  when  5  =  Pf/pP  ~  1.  Full  consideration  of  the  Basset 
term  is  especially  important,  but  difficult  from  a  numerical  point  of  view  owing  to  the  need 
of  storing  the  past  trajectories  of  a  large  number  of  particles.  Further  complications  arise 
from  considering  non-spherical  particles,  a  fact  which  induces  new  effects,  see  for  example 
[8]  for  experimental  results  on  the  dynamics  of  sinking  particles  with  plankton-like  shapes 
and  [5]  for  a  theoretical  study  of  the  dynamics  of  elongated  impurities. 
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(a)  <5  =  1.1 


(b)  5  =  0.9 


Figure  4:  Distribution  of  10,000  impurities  with  density  ratio  5  =  1.1  (left  panels)  and 
<5  =  0.9  (right  panels)  and  St  =  1  in  the  advecting  flow  defined  by  the  stream  function  ^  = 
2(cos  x+cos  y ).  The  individual  panels  refer  to  different  times,  namely  t  =  0, 0.5, 1,  2,  5, 10, 40 
and  100  non-dimensional  time  units. 


Figure  5:  (a)  Trajectory  of  an  individual  advected  impurity  with  density  ratio  6  =  0.9  and 
St  =  1  in  the  flow  defined  by  ^  =  2(cos  x  +  cos  y).  from  t  =  0  to  t  =  10000  non-dimensional 
time  units,  (b)  Diffusion  coefficients  along  x  and  y,  Dx  (stars)  and  Dy  (crosses),  for  an 
ensemble  of  heavy  impurities  as  a  function  of  e  =  1  —  5. 


Figure  6:  Streamfunction  and  the  corresponding  velocity  field  for  the  anisotropic  cellular 
flow  described  by  "0  =  (^4cosx  +  B  cos  y)  with  A  =  2  and  B  =  2.4. 


120 


Figure  7:  Distribution  of  10,000  impurities  of  density  ratio  5  =  0.9  and  St  =  1  in  the 
Eulerian  flow  defined  by  the  anisotropic  stream  function  ip  =  (Tlcosx  +  Bcosy).  A  =  2, 
B  =  2.4. 


(a)  (b) 

Figure  8:  (a)  Trajectory  of  an  individual  advected  impurity  with  density  ratio  <5  =  0.9  and 
St  =  1  in  the  Eulerian  flow  defined  by  the  anisotropic  stream  function  ip  =  (A  cos  x+B  cos  y ) 
with  ^4  =  2  and  B  =  2.4,  from  t  =  0  to  t  =  10000  natural  time  units,  (b)  Diffusion 
coefficients,  Dx  (stars)  and  Dy  (crosses),  for  heavy  impurities  as  a  function  of  e  =  1  —  5,  for 
A  =  2  and  B  =  2.2. 


Figure  9:  (a)  Fluid  parcel  trajectories  (thick  lines)  and  magnitude  of  the  Okubo- Weiss  pa¬ 
rameter  Q  (background  shading,  lighter  shading  correspond  to  positive  values  of  Q),  for 
the  stream  function  ^  =  A  cos  x  cosy,  (b)  Trajectory  of  a  neutrally  buoyant  impurity  with 
Stokes  number  St  =  0.2  (solid  line)  and  path  of  a  fluid  particle  starting  at  the  same  ini¬ 
tial  position  of  the  impurity  (dashed  line).  The  instability  in  regions  of  high  Q  allows  the 
impurity  to  wander  between  cells,  diverging  from  the  streamlines,  (c)  After  a  complicated 
excursion,  the  neutrally  buoyant  impurity  eventually  settles  in  a  zone  of  low  Q.  (d)  Dif¬ 
ference  between  the  velocity  of  the  impurity  and  the  fluid  velocity  at  the  instantaneous 
impurity  position. 
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Figure  10:  The  upper  panels  show  the  Poincare  sections  of  a  fluid  parcel  trajectory  (a)  and 
of  a  neutrally  buoyant  impurity  with  St  =  0.2  (b)  in  a  time-dependent  cellular  flow  with 
A  =  250,  B  =  0.3,  and  uj  =  1.0.  Panel  (c)  shows  the  time  series  of  the  difference  between 
the  velocity  of  the  impurity  and  the  fluid  velocity  at  the  impurity  position. 


a 


Figure  11:  Small  neutrally  buoyant  impurities  with  Stokes  number  St  =  0.2  tend  to  con¬ 
centrate  in  the  centers  of  vortices  in  a  two-dimensional  turbulent  flow.  Panel  (a)  shows  the 
distribution  of  the  Okubo- Weiss  parameter  Q  at  time  t  =  1  non-dimensional  units  (lighter 
shading  corresponds  to  larger  Q,  that  is,  to  hyperbolic  regions).  Panel  (b)  shows  the  dis¬ 
tribution  of  neutrally  buoyant  impurities  which  were  initially  uniformly  distributed  in  the 
flow  domain. 
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Figure  12:  Distribution  at  t  =  6  non-dimensional  time  units  of  8000  heavy  impurities 
advected  by  forced,  statistically  stationary  two-dimensional  turbulence  in  a  cyclonically 
rotating  reference  frame.  Small  solid  triangles  indicate  the  positions  of  the  advected  particles 
and  thin  solid  curves  indicate  vorticity  isolines.  The  impurities  concentrate  in  the  cores  of 
anticyclonic  vortices.  Parameter  values  are  <5  =  0.8,  St  =  0.01,  and  Ro  =  1/320. 


Figure  13:  Panel  (a)  shows  the  vorticity  field  at  t  =  15  non-dimensional  time  units  for  our 
idealization  of  the  flow  in  a  protosolar  nebula.  Panel  (b)  shows  the  distribution  of  heavy 
impurities  at  the  same  time. 
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Plankton  sinking  and  turbulence 

Lecture  given  on  1  July  2010  by  A.  Provenzale 


1  Introduction 

In  still  fluid,  heavy  particles  tend  to  fall  and  reach,  after  some  time,  a  terminal  velocity 
determined  by  the  balance  between  the  pull  of  gravity  and  the  viscosity  of  the  fluid.  Often, 
plankton  cells  are  slightly  heavier  than  sea  water,  and  as  a  consequence  they  tend  to  sink 
in  the  water  column.  This  sinking  could  then  take  them  below  the  euphotic  layer,  arresting 
photosynthesis. 

To  overcome  sinking,  plankton  have  devised  various  types  of  strategies.  One  is  swim¬ 
ming,  that  is,  an  active,  energy-consuming  self-propulsion  than  can  oppose  gravitational 
sinking.  Other  phytoplankters  can  modify  their  buoyancy,  becoming,  at  least  for  some 
time,  positively  buoyant.  Others  can  exploit  the  turbulence  in  the  mixed  layer,  using  it  to 
stay  suspended  for  longer  times.  This  latter  option  was  discussed  at  length  by  Margalef  [5], 
who  associated  the  level  of  turbulence  in  the  surface  layer  with  the  type  of  organisms  that 
are  environmentally  favoured.  The  terminal  velocity  of  spherical  impurities  with  the  same 
density  grows  with  the  squared  radius  of  the  falling  impurity.  Thus,  larger  plankton  cells 
tend  to  sink  faster  than  smaller  ones,  provided  they  can  be  approximated  by  a  somewhat 
spherical  shape.  When  the  turbulence  in  the  euphotic  layer  is  low,  Margalef  argued,  smaller, 
possibly  swimming  phytoplankters  such  as  flagellates  are  favoured  as  the  larger  species  sink 
more  rapidly  out  of  the  euphotic  zone.  On  the  other  hand,  when  the  level  of  turbulence 
is  high,  larger  plankton  can  benefit  from  turbulent  suspension  and  become  favoured  with 
respect  to  the  small  ones  (which  feel  diffusion  limitation  of  nutrients  more  severely). 

One  potential  problem  with  this  view  is  that  it  is  not  clear  whether  turbulence  really 
favours  a  prolonged  suspension  of  heavy  impurities.  For  example,  the  results  of  an  experi¬ 
ment  by  Ruiz  et  al  [8]  apparently  indicate  that  turbulence  make  heavy  particles  sink  faster 
than  in  still  fluid.  However,  other  experiments  considering  different  flow  configurations  gave 
the  opposite  result,  see  for  example  [1].  A  careful  study  of  this  issue  is  thus  required. 


2  Suspension  of  heavy  impurities 


The  non-dimensional  equation  of  motion  for  a  heavy  impurity  under  the  influence  of  gravity 
is  written  as 


(TV 

dt 


St 


(V-u  +  z  W) 


(1) 
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